Magnetic Impurity in Bernal Stacked Bilayer Graphene 
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We investigate a magnetic impurity in Bernal stacked bilayer graphene by a non-perturbative 
numerical exact approach. In the two cases we study, impurity is placed on the top of two different 
sublattices (A and B) in bilayer graphene. We find that similar to the monolayer case, magnet 
moment of the impurity could still be tuned in a wide range through changing the chemical potential. 
However, the property of the impurity depends strongly on its location due to the broken symmetry 
between sublattices A and B caused by the Bernal stacking. This difference becomes more apparent 
with the increase in the hybridization and decrease in the on-site Coulomb repulsion. Additionally, 
we calculate the impurity spectral densities and the correlation functions between the impurity and 
the conduction-band electrons. All the computational results show the same spatial dependence on 
the location of the impurity. 
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Following the fabrication of monolayer graphene 
(MLG)fi bilayer graphene(BLG) has attracted inten- 
sive attention due to its unique electrical flexibility^— 
and unusual physical properties including the unconven- 
tional quantum Hall effects^— At the same time, BLG 
is regarded as a good candidate in spintronics because 
compared with MLG, BLG shows longer spin-relaxation 
times at room temperature ; 1 1 ' 1 2 which is crucial for spin- 
tronic devices In particular, transition metal atom is 
usually used as a spin provider in spintronic devices^ 
and this motivates us to study the properties of mag- 
netic adatom in the BLG system. 

Because of the vanishing density of states (DOS) at the 
Dirac point, when we dope MLG with magnetic adatoms, 
their local moments could be conserved well at finite 
temperature^ 5 - The behavior of magnetic impurity in 
MLG has been considerably studied^— and Kondo ef- 
fects are difficult to observe in this system. Significantly 
different from the MLG, the two sublattices in Bernal 
stacked 2 ^ BLG are no longer equivalent. The DOS near 
the Dirac point would be largely modified 2 ^ because of 
the degeneracy lifted by the inter-layer hoppings and the 
local density of states (LDOS) have distinguishing values 
on the two sublattices which plays a more important role. 
This spatial inhomogeneity in LDOS could be directly 
detected experimentally with the scanning tunneling mi- 
croscopy (STM)^ In BLG with Bernal stacking, it has 
been reported that whether the missing carbon atoms 
(vacancy) are generated on sublattice A or B greatly al- 
ters the defect-induced magnetism and the localization of 
the zero mode in vacancy statesj 22 i 23 In addition, the use 
of STM tip provides the possibility to control the position 
of adatoms with atomic precision on the two-dimensional 
open surface i 24 ' 25 

In Fig. [TJ we show such asymmetry: the carbon atoms 
on sublattice A belonging to layer 1 lies directly on the 
top of those belonging to layer 2 so that the number of 
inter-layer hopping is one, while the carbon atoms on 
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FIG. 1: (Color online), (a) Lattice structure and various 
hopping energies of Bernal stacked BLG . (b) Top view of 
BLG with the black(dotted) line forming the top layer and 
the orange(dashed) line forming the bottom layer, ai and a2 
are the surface unit vectors. 



sublattice B in one layer is just fixed at the center of 
hexagon in the other layer so that it has three inter-layer 
hoppings. The direct consequence is that the LDOS in 
Bernal stacked BLG has spatial inhomogeneity. In the 
vicinity of Dirac point, the LDOS on sublattice A is much 
smaller than that on sublattice B. We find that the local 
moment of the magnetic impurity on top of sublattice 
A could develop better than that on sublattice B. This 
difference is not only reflected by the physical quanti- 
ties of impurity itself but also by the spatial correlation 
functions between the impurity and the carbon electrons. 

In this paper, we report a quantum Monte Carlo 
(QMC) study of a magnetic impurity placed on the top 
of the sublattices A and B in the Bernal stacked BLG 
and compare the results with their monolayer counter- 
part. The QMC method we used deals with infinite sea 
of conduction electrons and many-body effect without 
any approximations, so the results we present are essen- 
tially exact. The paper is organized as follows. In Sec. 
II, we give a brief introduction to our theoretical model, 
the Anderson impurity model in Bernal stacked BLG, as 
well as our numerical methods. We also compute LDOS 



2 



for sublatticcs A and B in BLG without impurity and 
compare them with the monolayer case. The purpose is 
to see the symmetry breaking generated by the Bcrnal 
stacking in bilayer system. In Sec. Ill, we present our 
results obtained from the quantum Monte Carlo (QMC) 
simulation. Firstly, we show the basic thermodynamic 
quantities on the impurity site, including occupancy, dou- 
ble occupancy, local magnet moments and spin suscep- 
tibility varying the chemical potential and the tempera- 
ture. Secondly, using the maximum entropy method with 
particle-hole symmetry, we do analytical continuation for 
the Green's function obtained from QMC to extract im- 
purity's spectral densities. Finally, we study the charge- 
charge and spin-spin correlation functions between the 
impurity and the carbon sites. All the calculations were 
worked out on both sublattices A and B and the results 
in bilayer case are compared with those in MLG. In Sec. 
IV, we summarize our results. 



II. MODEL AND METHODS 
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FIG. 2: (Color online). LDOS in the vicinity of Dirac point, 
from top to down is for sublattice B in BLG, MLG and sub- 
lattice A in BLG, respectively. The inset shows the detail 
of LDOS in the vicinity of Dirac energy. We set ti = 0.2t, 
t-i = O.li in BLG. 



Our starting point is the single-impurity Anderson 
model which has an impurity orbit of energy Ed and 
on-site Coulomb repulsion U ^ The impurity orbit hy- 
bridizes with a conduction band with strength V. The 
total Hamiltonian is 

H = Hq + Hi + Hi ■ 

Hq is a tight-binding Hamiltonian. For BLG it is 



H = -t a lma b mja + H.C. 

<i,j> ,m,(T 



(i) 



<hj> 

where a m i a (b m i a ) annihilates an electron with spin a at 
the site R m i a (Rmib) on sublattice A(B) of graphene's 
hexagonal structure in the m-th layer. The lattice struc- 
ture of BLG is shown in Fig. [TJ Intra-layer hopping t is 
the nearest neighbor hopping integral between sublatticcs 
A and B in the same plane and t is about 2.8eV~l Here 
t is used as the energy unit in our calculation. For the 
inter-layer hopping, we use Slonczewski-Weiss-McClure 
parametrization ; 28 ! 29 i.e., t\ « OAeV and £3 « 0.3eV and 
the configurations of these two couplings can be seen in 
Fig. (Ua). fi is chemical potential and is equal to zero in 
pure graphene. 

Hi is the impurity Hamiltonian 



Hi =^2(e d - v)d\d a + t/4vK' 



(2) 



where d a annihilates an electron with spin a at the im- 
purity orbit. 

Finally, H2 describes the hybridization between the im- 
purity adatom and a carbon atom on the BLG. If we 
place the impurity adatom on the top of a site in the 
first layer, then the Hamiltonian is written as 



if 2 =^(c 10j 4 + H.c). 



(3) 



If the adatom is added on the top of sublattice A, 
cioa = a W(T , otherwise ci a = hoa- 

Our principal computational tool is the Hirsch-Fye 
QMC algorithmp22. Basic thermodynamic properties of 
the impurity site are directly calculated by this algorithm 
and the spectral density is extracted by the method of 
Bayesian statistical inference. The Hirsch-Fye algorithm 
naturally returns the imaginary-time Green's function 
Gd(T) = J2a Gdo (t) of the impurity. With this Green's 
function, we can easily compute basic thermodynamic 
quantities of the impurity orbit, such as the expectation 
values of the total charge 



n d = (ndf + ridi), 
the local moment squared 

m d = ((«dt " n di) 2 ), 

and the double occupancy 

rid^ridi = (nd^ridi)- 



(4) 
(5) 
(G) 



According to the fact that the impurity charge can cither 
be zero or one, we note that 



m 2 d = n d -2n d ,n di . 



(7) 
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A non-zero value of m 2 d indicates the formation of a mo- 
ment on the adatom orbit. The closer this value is to 
one, the more fully developed is the moment. We also 
calculate the static impurity spin susceptibility 

X = / dT(m d (T)m d (0)), (8) 
Jo 

where j3 = T _1 , m d [r) = e TH m d (Q) e - TH . 

Using imaginary-time Green's function obtained from 
the QMC method, we can calculate the spectral density 
A(w) = A a (uj) by numerically solving^! 

oo 
— oo 

The detailed procedure of Bayesian inference method is 
presented in Ref. [3ll This Bayesian inference procedure 
is also called the maximum entropy method. 

We also use an extension of the Hirsch-Fye algorithm^ 
to compute the charge-charge correlation function 

Ci = (nam) - (nd)(rii), (10) 

and the spin-spin correlation function 

Si = (mdrrii). 

Using the standard particle-hole transformation 
on one of the sublattice, we can prove that 
Gd(r,fi) — Gd{—T,—(i) and consequently that 
A a (/i, oj) = A a (—(j,, — oj) when e d = —U/2M These two 
results in turn imply the symmetries with respect to the 
sign of n in various thermodynamic quantities of interest. 
These symmetries also mean that without loss of gener- 
ality, we can restrict our attention to the behavior of the 
system when /i < 0. 

The DOS of MLG and BLG are distinguishing in the 
vicinity of the Dirac point because of their different dis- 
persion relations^. In general, the degeneracy near the 
Dirac energy are lifted due to the inter-layer hoppings *i 
and <3, so BLG shows larger DOS than MLG near the 
Dirac energy^ To gain some primary insights of the spa- 
tial inhomogeneity in BLG, we first calculate the LDOS 
without the magnetic impurity (V = 0). The results are 
shown in Fig. [2] as well as the LDOS for monolayer case. 
In MLG, every site is identical so that the DOS is equal to 
the local one while in BLG the asymmetry between two 
sublattices occurs. In particular, as is shown in Fig. [5J it 
is clear that the LDOS on sublattice B has a larger value 
than that on sublattice A and the value of monolayer is 
in between. As a result, we can expect that the differ- 
ence in LDOS for two sublattices will lead to position- 
dependent effective hybridizations between impurity and 
carbon atoms, and the features of impurity have spatial 
inhomogeneity which can be seen from QMC results in 
the following sections. 
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FIG. 3: (Color online), (a) rid as a function of fi. (b) 
ndupiT-ddown as a function of fi. (c) m\ as a function of fj,. 
(d) x as a function of [i. In all plates, U = 1.6*, Ed = —U/2, 
fj = 1/T = 40*~\ V = 1.0*. MLG: impurity added on top of 
MLG; A: impurity located on top of sublattice A in BLG; B: 
impurity located on top sublattice B in BLG. 
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FIG. 4: (Color online), (a) The spin susceptibility x as a 
function of T -1 for case A . (b) x as a function of T -1 for 
case B . U = 1.6*, e d = -U/2, V = 1.0*. 



III. RESULTS 
A. Basic Thermodynamic Quantities 

In Fig. [3] wc present the results of several thermody- 
namic quantities on impurity site as a function of fj, . 
The error bars are smaller than the points except where 
shown. In the following sections and figures, the case 
MLG, A and B is referred to the impurity located on 
top of MLG, sublattices A and B in BLG, respectively. 
Since e d = —U/2 and 11 = 0, the system has particle- 
hole symmetry. This symmetry fixes the total charge 
rid = {n d f + Tidi) exactly at one in Fig. [3](a) . Seeing from 
Fig. [3J it is clear that the four quantities decrease when 
/i moves below the Dirac point. In the values of charge, 
there is an order: case A > case MLG > case B . The 
values of the local moment squared m d and the spin sus- 
ceptibility x a l so show the same order, while the double 
occupancy shows the opposite order to them. This re- 
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FIG. 5: (Color online). m d as a function of fx for different 
V, from top to down is V = 0.75t, l.Ot, 1.25t. In all cases, 
U = 1.6t, e d = -17/2, p = l/T = 40t _1 . 
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FIG. 6: (Color online). m d as a function of fj, at U = 0.8t for 
different V, from top to down is V = 0.75t, l.Ot, 1.25t. In all 
cases, e d = -(7/2, /3 = 1/T = 40t~ 1 . 



suits from the relation 2nd U p?T-ddown = n<i — fn% which is 
satisfied automatically in our simulation although these 
three quantities are calculated independently. In Fig. |H 
we study the spin susceptibility x(^) versus temperature 
T in cases A and B . In both cases, as /i moves below 
the Dirac point and T is lowered, we see that \ crosses 
over from a Curie- Weiss behavior to the T-indcpendcnt 
behavior as a screened local moment. In particular, near 
the Dirac point, we see that the local moment on the 
sublattice A is developed better that on the sublatticc 
B as the temperature is lowered. In principle, spin sus- 
ceptibility defined in Eq. (|8) depends not only on the 
local moment itself but also on the spin correlation with 
conduction-band electrons^ As a result, when the LDOS 
of conduction electrons is higher, the spin correlation is 
larger, and \ nas smaller value and this is consistent with 
the results of LDOS in Fig. [2] Later we will show spin cor- 
relation directly, and this point can be seen more clearly. 

In order to see how hybridization V could affect the 
local moment, we do the calculation with different values 
of V and the results are shown in Fig. [5] We can see 
that as V grows from 0.75t to 1.25t, the local moment 
of the impurity site decreases noticeably for all the three 
cases we examine. As V varies, the order for the three 
cases shown in Fig. [3] does not change. The Coulomb 
interaction U for transition-metal atom in carbon-based 
materials can be varied from 2-5 eVj2£r— so we also study 
the effects of different on-site Coulomb repulsion. Shown 
in Fig. [S] are the results for U = 0.8t with other parame- 
ters be the same as those in Fig. [SJ We can see that the 
general behavior of the local moment remains the same 
with Fig. [5] as we vary \i and V, and the order for the 
cases MLG, A and B persists. Comparing Fig. [5] with 
Fig. [13 we see that the distinction for the cases A and B 
are more obvious with a smaller U. 
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FIG. 7: (Color online). A(cj) with respect to u) for [a)V = 
0.75i, U = 1M, (b)V = l.Ot, U = 1.6* and (c)V = l.Ot, 
U = 2 At. Here e d = -(7/2, f3 = 1/T = 40t _1 in all cases. 



B. Spectral Densities 

In Fig.[7]are the spectral densities A{u) at 1/T = 40t _1 
with various U and V. We fix fi = and = —(7/2, 
so there is particle-hole symmetry and we have A(— cj) = 
A(uj). It is well known that for the Hartrcc-Fock solu- 
tion of an Anderson impurity in a normal metal^ the 
locations of the peaks arc independent of V and the sep- 
aration of two peaks is about U. However, here we see 
that the separations of two side peaks in A(ui) are much 
smaller than Coulomb repulsion U, this point is con- 
sistent with the results from numerical renormalization 
group studyJ^ When V increases, the A(u) peaks move 
toward the Dirac point. These results are in agreement 



with those in the previous study at higher temperature^ 
in the case MLG , which is absolutely originated from 
the non-constant DOS in graphene. We can also find 
that near the Dirac point, for the three cases studied, 
the spectral densities have the same order as that of the 
LDOS shown in Fig. [2] near the Dirac point. This is 
because when V is induced, the eigenstates in host sys- 
tem can greatly hybridize with impurity orbit, and as 

V increases, the LDOS on two sublattices in graphene 
influences the impurity orbit more strongly. 

In Fig. [7£a) where V = 0.75i and U = 1M, we can see 
that the differences in spectral densities in the three cases 
are relatively small near the Dirac point. In Fig. (7{b), as 

V is increased from 0.75t to l.Ot, the peaks in A(oj) move 
towards the Dirac point and thus the A(ui) for case B 
increase dramatically while no such changes are induced 
in the other two cases. 

We also study the A(ui) with a larger value of Coulomb 
repulsion U = 2 At in Fig. 0c) and V = l.Ot with 
particle-hole symmetry. It is clear that the separations 
of two peaks in A{ui) increase but are still much smaller 
than U. Near the Dirac point, the spectral densities for 
the three cases have the same order as that in Fig. [TJa)- 
(b). 



C. Correlation Functions 

In Fig. [51 we present the results of charge-charge cor- 
relation Ci and spin-spin correlation Si between impu- 
rity and conduction- band electrons. We set V = l.Ot, 
U = 0.8t, £d = —U/2. The chemical potential is fixed 
at zero point, so the system has particle-hole symmetry. 
In every subfigure, the impurity is located on the top 
of the site i = 0, so the locations with even index are 
sites on the same sublattice as the site where impurity 
is added and those with an odd index are sites on the 
opposite sublattice. As shown in Fig. [5J Ci and Si for all 
the three cases are relatively short-ranged such that the 
magnitudes decay rapidly with respect to the location i. 
Ci lacks oscillations since at fj, = the total system is 
half filled and the charge exchange between two sites is 
greatly suppressed. If we look at the correlations in de- 
tails from the insets (ai) and (a 2 ) in Fig. [5J we find that 
compared to the case B, the on-site correlation |C 1 of the 
case A is weaker while the nearest-neighbor correlation 
Ci is stronger. 

The plots of Si show that at /i = 0, the impurity spin 
is antiferromagnetically correlated with the conduction 
electron spins on the same sublattice, and fcrromagncti- 
cally correlated with those on the opposite sublattice. To 
compare Si in the cases A and B , we also present the 
insets (bi) and (62), and we see that So is stronger for 
the B case while ferromagnetic correlation Si is stronger 
for the case A , which is consistent with Ci. Due to the 
bipartite nature of the lattice and the localization of im- 
purity spin, the system would have weak spin fluctuation 
around the adatom at half filling. 




FIG. 8: (Color online). (a)The charge-charge correlation d 
and (b) the spin-spin correlation Si versus site i. i is along a 
zigzag direction in layer 1 of BLG, where impurity is located 
on top of site i — 0. Insets (ai) and (02): details of Co and Ci; 
insets (bi) and (62): details of So and Si. In inset (63), the 
sites along red(dashed) and blue(dotted) zigzag lines indicate 
the carbon sites we consider for case A and B, respectively. 
The two red circled sites are the sites i = 0, where impurity 
is located for cases A and B. Here we use V — l.Ot, U — 0.8t 
and E d = -U/2, /3 = 1/T = 40£ _1 in all the cases. 
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FIG. 9: (Color online). (a)The charge-charge correlation Co 
and (b)the spin-spin correlation So versus /i. Here we use 
V = l.Ot, U = 0.8t and e d = -U/2, /3 = 1/T = 40f"\ 



We also focus on the behavior of on-site correlation 
functions Co and Sq with (i moving below the zero point 
in Fig. |5l The particle-hole symmetry is broken and the 
filling is shifted from one. Here V = l.Ot, U = 0.8t , 
e c i = —U/2 and j3 = 40£ -1 . As fi is tuned below the 
Dirac point, the amplitude of on-site charge correlations 
Co in Fig. Uta) increase because the occupancy of the 
impurity orbit and conduction band shift from half filling 
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so their charge exchange enhances. The differences of Co 
for the three cases are the most obvious at [i = and 
as fi is lowered, the differences become smaller and at 
/x « — 0.2t, three curves touch each other. We can see 
the same behavior of So in Fig.^b) as that in Cq. 

We can correlate the Si with the spin susceptibility in 
Fig. GJd) and Fig. At /i = with particle-hole sym- 
metry, all of the three cases we study have half filling, so 
there exists well-defined local moment on impurity site. 
The main differences for spin susceptibilities among them 
originate form the screening of conduction-band electrons 
(in fact T\ is screened moment), which is reflected by the 
spin-spin correlations. Furthermore, Si is directly de- 
pended on the LDOS of conduction band in Fig[2j so in 
Bernal stacked BLG, the spatial inhomogeneity for spin 
susceptibility can be understood. 

IV. CONCLUSIONS 

In summary, we have studied a magnetic impurity 
placed on the top of two nonequivalent sublattices in 
Bernal stacked BLG with Slonczewski-Weiss-McClure 
parameterization. The results obtained from the quan- 
tum Monte Carlo method are essentially exact, in the 
sense that we start with an infinite sea of conduction 
electrons and use no approximations to deal with many- 



body problem in our simulations. The LDOS on the 
two sublattices show spatial inhomogeneity in BLG. As 
a result, when we put magnetic impurity on the top of 
the two sublattices, such spatial inhomogeneity greatly 
influences the magnetic property of the impurity. It 
is interesting that this inhomogeneity was also seen in 
vacancy-induced magnetism in BLG experimentally^ In 
general, we found that the local moment on the sublat- 
tice A is conserved better than that on sublattice B and 
this difference becomes more apparent as hybridization 
V increases and Coulomb repulsion U decreases. Other 
physical quantities mostly have the same feature. The 
STM could be used to measure the spectral densities 
and the charge-charge correlation functions, and a spin- 
polarized STM could be used to measure the spin-spin 
correlations ; 39-41 so our studies are well connected to ex- 
periments. 
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